﻿#include <algorithm>
#include "../common/common.h"

const int N = 33 * 1024;
const int threadsPerBlock = 256;

__global__ void dot( float *a, float *b, float *c ) {
	__shared__ float cache[threadsPerBlock];
	int tid = threadIdx.x + blockIdx.x * blockDim.x;
	int cacheIndex = threadIdx.x;
	float temp = 0;

	while (tid < N) {
		temp += a[tid] * b[tid];
		tid += blockDim.x * gridDim.x;
	}

	// ghi kết quả vào bộ nhớ chia sẻ
	cache[cacheIndex] = temp;

	// đồng bộ hóa các tiểu trình trong block hiện tại
	__syncthreads();

	// với bài toán reduction, threadsPerBlock phải là 
	// lũy thừa của 2 để đoạn mã này chạy đúng
	int step = blockDim.x/2;
	while (step != 0) {
		if (cacheIndex < step)
			cache[cacheIndex] += cache[cacheIndex + step];
		__syncthreads();
		step /= 2;
	}

	if (cacheIndex == 0)
		c[blockIdx.x] = cache[0];
}

const int blocksPerGrid =std::min( 32, (N+threadsPerBlock-1) / threadsPerBlock );

int main( void ) {
	float *a, *b, c, *partial_c;
	float *dev_a, *dev_b, *dev_partial_c;

	// cấp phát bộ nhớ trên CPU
	a = (float*)malloc( N*sizeof(float) );
	b = (float*)malloc( N*sizeof(float) );
	partial_c = (float*)malloc( blocksPerGrid*sizeof(float) );

	// cấp phát bộ nhớ trên GPU
	SAFE_CALL( cudaMalloc( (void**)&dev_a,
		N*sizeof(float) ) );
	SAFE_CALL( cudaMalloc( (void**)&dev_b,
		N*sizeof(float) ) );
	SAFE_CALL( cudaMalloc( (void**)&dev_partial_c,
		blocksPerGrid*sizeof(float) ) );
	
	// ghi giá trị cho dữ liệu trên CPU
	for (int i = 0; i < N; i++) {
		a[i] = i + 10;
		b[i] = i * i;
	}

	// chép mảng 'a' và 'b' từ CPU sang GPU
	SAFE_CALL( cudaMemcpy( dev_a, a, N*sizeof(float),
		cudaMemcpyHostToDevice ) );
	SAFE_CALL( cudaMemcpy( dev_b, b, N*sizeof(float),
		cudaMemcpyHostToDevice ) );
	
	dot<<<blocksPerGrid,threadsPerBlock>>>( dev_a, dev_b,
		dev_partial_c );

	// chép mảng 'c' từ GPU về CPU
	SAFE_CALL( cudaMemcpy( partial_c, dev_partial_c,
		blocksPerGrid*sizeof(float),
		cudaMemcpyDeviceToHost ) );

	// tính toán bước cuối cùng trên CPU
	c = 0;
	for (int i=0; i<blocksPerGrid; i++) {
		c += partial_c[i];
	}

#define sum_squares(x) (x*(3*x*x*x + 46*x*x + 63*x + 20)/12.0)

	printf( "Does GPU value %.6g = %.6g?\n", c,
		sum_squares( (double)(N - 1) ) );
	
	// giải phóng bộ nhớ GPU
	cudaFree( dev_a );
	cudaFree( dev_b );
	cudaFree( dev_partial_c );
	
	// giải phóng bộ nhớ CPU
	free( a );
	free( b );
	free( partial_c );
}